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1. Introduction 

Models of finitely connected disordered spin systems have been studied for some twenty 
years, following the initiating papers [U El El IH Ej • Especially due to the unexpectedly 
rich and varied range of multi-disciplinary applications of finite connectivity replica 
techniques which emerged subsequently, in e.g. spin-glass modeling [3131313, error 
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correcting codes (TUl HU C21 EE] , theoretical computer science [EJ [T51 EH U7j , recurrent 
neural networks fHl El 120]; and 'small- world' networks |21j . this field is presently 
enjoying a renewed interest and popularity. Until very recently, analysis was limited 
to the equilibrium properties of such models, but now attention has also turned to the 
dynamics of finitely connected spin systems [221 ESI EH ESj , using combinatorial and 
generating functional methods. In the domain of physical spin systems, research into 
finitely connected systems has usually been triggered by the desire to develop solvable 
spin-glass models which are closer to real finite-dimensional systems than the celebrated 
fully connected spin-glass model of [26J. As far as we are aware, however (and in contrast 
to the situation with fully connected disordered spin systems), all finitely connected 
and disordered spin systems analyzed theoretically so far involved scalar spin variables 
(either of the Ising type, the soft-spin type, or the spherical type). 

In the present paper we solve equilibrium models of finitely connected spin systems 
of unit-length vectorial spins, and with random pair-interactions between them which 
are of a chiral nature, defined by random orthogonal matrices which promote random 
relative spatial orientations between pairs of spins. The motivation behind our study is 
twofold. Firstly, we aim to expand the domain of solvable and solved finitely connected 
spin models, by also including those where the spins have a truly vectorial character, 
and where their interactions are of a (random) matrix type. Vectorial spins are not only 
more realistic from a fundamental physical point of view in any magnetic system, but 
are also of special relevance in the context of e.g. Josephson junctions [23 [2H1 121] • In 
these latter junctions it is essential for the spin-interactions to have a chiral character. 
The inclusion of bond- and field-disorder in chiral spin systems has so far only been 
studied analytically in a mean-field setting, where all spins are allowed to interact 
with each other (purely for mathematical convenience, see e.g. [30J and references 
therein). The present study can be regarded as a step away from the unrealistic full 
connectivity towards finite dimensional models in such vectorial systems. Our second 
motivation is of a technical nature. In contrast to models with discrete (e.g. Ising) spins, 
when choosing real- valued microscopic variables the replica-symmetric (RS) theory will 
involve an order parameter which is itself a functional, rather than a function. Solving 
the associated order parameter equations is therefore non-trivial, and leads to many 
numerical complications especially when the domain of values for the individual spins is 
not bounded (as is the case for e.g. soft spins). In our present model we have vectorial 
spin states; although continuous and therefore leading to a theory involving an order 
parameter functional, each spin state represents a point on a sphere and has therefore 
a compact domain. This is found to be a considerable mathematical and numerical 
advantage, and allows us to push our analysis and therefore also our understanding 
further. 

The structure of our paper is as follows. We first develop the general 
replica symmetric theory for finitely connected vectorial spin systems with random 
chiral interactions described by orthogonal matrices, for arbitrary dimensions of the 
microscopic sphere Sd-i which constrains the values of the individual spins. We then 
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apply our theory first to the case d = 2, where the spins reduce to XY ones, followed by 
application to d = 3 where we have classical Heisenberg spins. In both cases we calculate 
phase diagrams (with continuous transition lines calculated via functional moment 
expansions) as well as the values of macroscopic observables (using truncated population 
dynamics routines), for different choices for the ensemble of random orthogonal matrices. 
We complement our results with numerical simulations. The latter are found to agree 
well with our theoretical predictions, especially if the usual problems associated with 
the moment-truncation of population dynamics algorithms and the well-known finite 
size effects in disordered spin systems are taken into account. 

2. Definitions 

We study finitely connected and bond-disordered systems composed of iV interacting 
normalized vectorial soft spins <Ti e Sd-i, with Sa-i denoting the unit sphere in IR d . 
Thus for d — 2 our spins become XY spins, with each spin representing a point on a 
the unit circle, whereas for d = 3 they become classical Heisenberg spins, with each 
representing a point on the unit sphere, etc. Our systems are taken to be in thermal 
equilibrium, characterized by Hamiltonians of the following form, with the short-hand 
{cr} = ((T 1 ,...,(T N ): 

H({CT}) = -J£c^ • V ij( Tj + ^V(<T i ) (1) 

i<j i 

The independently distributed quenched random variables c^- G {0, 1} define the 
(random) connectivity of the system, i.e. they dictate which pairs of spins are allowed 
to interact. The real d x d matrices are assumed orthogonal and with determinant 
one, i.e. they represent rotations in III d , and are drawn randomly and independently 
from some random matrix ensemble characterized by a distribution P(U). We will 
assume as yet only that P(U) = P(U^). The single-site potentials V(cr) serve to break 
symmetries, and will also enable us to take an Ising limit later as a benchmark test 
(e.g. for V(cr) = 7(e • cr) 2 , with e denoting a fixed unit length vector in IR d , and 
where 7 — > 00). In this paper we will be concerned with the so-called finite connectivity 
regime, where 

P(cij) = ^6 Cij ,i + (1 - ^)S Cij ,o for all i < j (2) 

with c = O(N ). In this limit, each spin interacts on average only with a finite number c 
of other spins, even in the thermodynamic limit (similar to finite dimensional systems). 
In this sense, in spite of the absence of geometrical considerations, finitely connected 
spin models can be regarded as closer to physical reality than fully connected ones. 

In the remainder of this paper, averages over the random connectivity {%} and 
over the random ortoghonal matrices {U^} will be denoted by 7TT . We will use the 
standard replica identity logZ = lim n _>o ^ 1 log Z n , with Z n initially evaluated for 
integer n, to calculate for our system the asymptotic disorder-averaged free energy 
per spin / = — lim A r_ >00 (/3A r )~ 1 log Z. Indices will be used according to the standard 



Finitely connected vector spin systems with random matrix interactions 



conventions, with Greek ones labeling replicas (a = 1, . . . ,n) and Roman ones labeling 
spins (i — 1, . . . , N). 



3. Replica calculation of the disorder-averaged free energy per spin 

3.1. Derivation of replica saddle-point equations 

In evaluating the disorder-averaged free energy per spin with the replica method, 
1- 



lim 



log/ []J da i] e-? H (W) 



one finds that site factorization can be achieved upon isolating the usual site- 
averaged replica order parameter for finitely connected systems. After redefining 
{cr} = (cr 1 , . . . , cr n ), where now cr°,cr a 6 Sd-i, this order parameter takes the form 
P({cr}) = limjv^oo N" 1 J2i Ha S[cr a — erf]. The insertion of appropriate functional 6- 
distributions into (J3J) gives 

7= - lim lim -i- log [ TT dtrfe-fWV] TT ^ ^f^f 



lim lim — 

N-oonr+O /3Nn J JS d ^ ,„ 



log /. . . 



[lid 



x exp j Jl^ [| dU p( U )e' 3J S Q ^-u^ _ i] + (N°)\ 

, i 1 i f f tt dP({<7})dP({ CT })e iJVP « CT » /5 «°'» 
= — hm hm log / < — — — — — — —, 

v^oo n^o (3Nn B J \ 2vr /N 

x expi^cN J{d(Tdcr'} P({(r})P({tr'})[j "dU P(U)e /3J S a °V ucr '« - 1] 
x exp [N log /{d<x} e -»E {< r > ^'})IL^-'-^]-/»E„v(*r.)j 

In this expression we substitute P({<x}) — > {dcr}P({cr}); this subsequently enables us 
to take a continuum limit, whereby {dcr} — > {0}. Upon writing the resulting path 
integration measure as Il{cr}[dP({cr})dP({cr})/27r] = {dPdP}, and upon neglecting 
irrelevant constants, our expression for the disorder averaged free energy per spin is 
then seen to take a saddle-point form: 

7 = - lim lim-!- log [ {dPd p } e iN f{d(T} P({CT})P({<T}) 

J n^oo n^o BNn J 1 J 



X 



exp \ ^cN j{dada'} P({ct})P({ct'}) [ j dU P(U)e /3J ^. <r a u<x-_ j 



x expjiVlog J{da} e -^E ^)-^(W)| 

- lim-^-extr {P p } jij{d<r} P({<x})P({<x}) + log j{d<x} e -*E a v(«r«)-i*{<r}) 
+i C |{d<7d<r'} P({<t})P(K}) [|dU P(U)e^E^ a -uo- a '_ i]| ( 4 ) 
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The functional variation of (@J) with respect to the two order parameter functions P({cr}) 
and P({cr}) gives us the following two saddle-point equations, respectively: 

P ^ CT ^ ~ J {da'} e-f^E^i "")-^ "}) ^ 

P({a}) = ic j{da'} P{{*'}) [ J dU P(U)e /3J S^ vaa '- l] (6) 

Elimination of the conjugate order parameter function P({cr}) leads to a saddle-point 
equation for P({er}) only, and an associated expression for the free energy /: 

e c J{d<T'} P({CT'})[/dU P^e^^'^^-l^^Vi^) 

PiW}) = fiT x- o-°"i,n"»" ( 7 ) 

/{do-'} e c / {d<T " } p({<T " })[ / dU p(u)e Q -i]-/?£ a v(o-') 

1 = H(L{l C ^(W)^(K})[/dU P(U)e^ Q ^-u<T«'_ 1 

-log|{do-} e c ^ do "} p ({°"})[/du P(u)/ J £ Q CTa uCTtt V/3E«n^)J ( 8 ) 

c?.^. Replica symmetric theory 

We now make the canonical RS ansatz for continuous spins in our saddle-point 
equations! . We assume there to be a complete family of distributions P[<r|/x] on S^-i, 
parametrized by a countable set of real-valued parameters /x = (//q, 1^2, ■ ■ •)> sucn 
that 

P RS (o- 1 ,...,o-")= /"d/i «;(/*) JJP[o- a |/i] (9) 

J a 

with a normalized density w(i-i). A representation- independent but mathematically 
equivalent formulation of the RS ansatz follows upon defining the functional measure 

W[{P}]= fd^w(ti) J] 6[P(tr)-P[<r\ri\ (10) 

For non-degenerate parametrizations, i.e. for those such that every function P(cr) 
corresponds to a unique choice /j,({P}) of parameters, with P[(t\h({P})] = P(o-) for 
all cr G Sd-i, we can invert relation (fTU|) and write 

w{n) = /{dP} W[{P}} % - M ({P})] (11) 

In terms of the functional measure ^[{P}] our RS ansatz Q takes an elegant and 
representation-free form: 

P RS (o-\...,o-)= /{dP}^[{P}] U P M (12) 

| This ansatz reflects the complicating fact that, in the case of continuous spins, the RS order parameter 
function depends on the replicated spin variables not only via the sum J2 a a a ( as wou ld have been the 
case for Ising spins), but rather on all possible sums of the form J2 a a ai f° r am/ K > 1. 
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The physical interpretation of our subsequent observables and results in terms of the 
original disordered iV-spin system will follow from the identity 

/{dP}^[{p}]n[/d^p( CT )/ a ( CT )l = hm ixinw) ( i3 ) 

We insert into our general saddle-point equation (|2J) the RS ansatz (fEZj) . and introduce 
the convention l\k=i a k = 1 for any series {a k }. This leads to the following identity 
(with a constant C n which will in due course be determined by normalization): 

C n [{dP} W[{P}]UP{<r a ) 

J a 

cf{dP}w[{P}] JdU P(u)n Q [/dcr'p(cr')e' 3 ^- uCr 'j- c -/3X; ci V'(cr ct ) 



E/ _c / nK^}^[{p,}]du fc p(u fc )] 

e>o l - J k=i 

a { k=l J J 

{dP}n^K)E^- c f{[[{dp k }w[{p k }]dv k p(v k )] z-[{p 

■ " J fe=l 



x n * 

(TeS d 



£>0 

P(<r)- 



Z[{P u ... i P i }] 



where 



Z[{Pu ■ ■ ■ , Pe}} = /dcr e-^W J[ / d ff 'P fe (tT') e ^' u ^' 

fc=i 



^}] 
(14) 

(15) 



In the replica limit n —>■ 0, both the term Z n [{Pi, . . . , P&}] and the constant C n reduce 
to unity, and our RS order parameter equation acquires the transparent form 

w[{p}\ = EtW ]j[{dp k }w[{p k }]dv k p(v k )] 



x n s 



k<l 

P(a) 



e -M<r) ni =1 /d<x'P A V)e /3J<T ' Ufc<T ' 



/d<r"e-^(°"") nLi Jd<x'P fc (cr')e/ 3Jcr ''- u * cr '_ ' U>) 

The replica symmetric order parameter for finitely and randomly connected disordered 
systems with continuous degrees of freedom is thus seen to be a functional ^[{P}] 
on the space of probability densities. For any specific parametrization P[<x|/x] of this 
space, the order parameter equation (jTfij) becomes an equation for the distribution w(fj,) 
of parameters (of which there must generally be an infinite number): 



w 



>) = E|^ C /n[ d M^WdU fc P(U fc )] J{dP}8[» - M ({P})] 

e -/3V((T) ^ / d0 -'P[ -'|/x fc ]e /3J<r - Ufc0 " 



k<t 



x n * 



P(tr) 



;i7) 



We note that multiplication of P[<r|/Lt] by a constant will not affect ([17)1 . so that in our 
parametrizations P[<r|/x] we need not impose normalization explicitly. Application of 
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the above manipulations to formula (JSJ) leads us in a similar manner to the following 
expression for the RS free energy: 

c 



/rs = ^ y{dPidP 2 } W[{Pi}]W[{P 2 }] ydUP(U)log 
3 E J^ C I LI [{dP fc }^[{P fe }]dU fc P(U fc )] 



|d -d -'P 1 ( -)P 2 ( - / )e /3JO '- UO " 



/3 



x 



logj^do- e-P y((J) n|d < T / P fc ( < T / )e /3J<T Ufe °"| (18) 



5.5. T7ie Ising limit 

As a simple consistency test of our theory, we now consider the limit 7 — > 00 of our 
population dynamics equation for the special choice V(<x) = 7(6. cr) 2 where e 

denotes an arbitrary unit-length vector in IR^. Via a saddle point argument, we observe 
that the limit 7 —>■ 00 restricts our spins to take one of two possible values, viz. cr = ±e. 
Hence, our order-parameter function P(er) collapses to a sum of two delta peaks, which 
(with a modest amount of foresight) we choose to parametrize by 



P(tr) 



-6(cr - e) + 



-j3h 



-5(cr + e) 



2cosh03/ip~ ~' ' 2cosh(/3/ip~ ' ~' ^ 

With this representation of P(cr), our order parameter functional M^[{P}] reduces to 
a distribution W(h) of 'effective' fields h. Carrying out the integrations within the 
5-functional in (j!6)l then results in 



W{h) 



E 

£>Q 



—c £ 

e c 

~1T 



]J[dh k W(h k )dV k P(U k ) 

k<£ 



(20) 



x n s 

CTe{-e,e} 



e /3h(CT.e) 

2 cosh(/3/i) 



nLi 


' e /J/ Ife + / 3jCT.U fe e e -/3h fe -0jCT.U fc e" 




2cosh( / 3/i i: ) 1 2cosh(/3fe fe ) 


Scr'e{-e,e} ]lfe=l 


e f3h k +/3j(T'.U k e e -f3h k -0j(T' .U k e 


2coshQ3h fe ) ' 2cosh(/3/i fc ) 



The two 5-functions in this expression now effectively give us an update relation for h. 
This leads to a population dynamics equation, equivalent to that found in e.g. [2HH|: 



— c £ 

W ( h ) = Y,^JT [ Hidh k W(h k )dU k P(U k )] 



£>0 



x 5 



h — — arctanh[tanh(/3/i) tanh(/3Je. U^e)] 



fc=i 



(21) 



Choosing, for instance, the random orthogonal matrix distribution P(U) to be of the 
simple form P(U) = aP(l) + (1 — a)P(— 1)§, with a e [0,1], leads directly to the 
familiar functional order parameter equations of the ± J Ising spin-glass, with exchange 
interactions distributed according to P(J') = a5[J' — J] + (1 — a)5[J' + J]. 

§ Here we take the freedom to consider the orthogonal group 0(3) instead of SO(3), otherwise the 
matrix 1 would not be in our ensemble (defined as the orthogonal matrices with determinant one). 
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4. General theory for d = 2 

4-1- RS saddle-point equations and free energy 

For d = 2 the set S^-i reduces to the unit circle, and our spins become XY-spins, i.e. 
er = (cos 0, sin 0) with G [0, 27r]. We may therefore also write V(cr) = V{(j)), and find 
the random rotation matrices U being simply characterized by a single angle u G [0, 2tt] 
and an associated symmetric distribution P{oj) = P(—u). As a result our equations 
simplify considerably. The order parameter equation (fTKj) reduces to 

w[{p}\ = E^ C [mww[{p k }]du k p(u k )} 

e -/3V(<t>) Yf k=1 fd(f)'P k ((j)')eP Jcos ^'- w ^ 



k<l 



p(0)- 



Jd0"e-^ y ^") n fc= i Jd^'Pk^ei 3 ^ ^"-*'-^ 
In the absence of single-site potentials, i.e. for V((p) = 0, we get 



(22) 



^>o 



p(0)- 



nLi Jd0'P fc (0 , )e /3Jcos{0 " </, '"" fc) 



For T = this reduces even further to 

nti^(0-^) 



-<f>'-ui k ) 



(23) 



£>0 



p(0) 



(24) 



Similarly we find that for V(cr) = and d = 2 the RS free energy per spin (j!8j) becomes 



./ 



RS 



c 

2^ 



_ 1 



{dpdP 2 } W[{Pi}]W[{P 2 }] yd^PM log 
n[{dP fe }M{P fe }]du; fe P(u; fe )] 

d0 / d0 / P fc (0 / )e^ JcoB( ^-*'- Wh ) 



d0d0'P 1 (0)P 2 (0 / )e /3JcosW '- < ^ ) 



fe=i 
x log 



fc=l • 



(25) 



As expected, the paramagnetic state W[{P}] = n</,e[o,27r] ^[P(4") ~ (27r) -1 ] is a solution 
of at any temperature. In this state one finds, using (fT3j). that (5[0 — 0j]) = (27r) _1 
for all i and all 0. 

Continuous bifurcations away from the paramagnetic state can be identified via a so- 
called Guzai (or functional moment) expansion. We transform P(0) -> (27T)" 1 + A(0), 
with W[{P}] -> W[{A}] and with W[{A}] = as soon as / O 2 M0 A(0) ^ (since 
P(0) must remain normalized). We may now expand our equations in powers of the 
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functional moments J{dA}W / [{A}]A(0 1 ) . . . A(0 r ) for r = 1,2. In doing so we will 
repeatedly encounter the modified Bessel functions 

U*)= r^cos(n0)e zcos W (26) 

Close to a continuous phase transition we assume there to be a small parameter e 
measuring the bifurcation, such that / {dA}W / [{A}] A(0i) . . . A(0 r ) = 0(e T ). 

4-2. Paramagnetic to ferromagnetic and Kosterlitz-Thouless type transitions 
In lowest order e 1 we have 

/{dA}#[{A}]A(0) = ^^4- c /nK d M^[{A fe }]d^PK)] 



£>0 *" fc<£ 



X 



1 + ^*=1 /o(/3J) + • • • 

14 _V^ rd^" /d^A fc (^) e g Jc "(»"-* , -'fc) 

1 T <Lfc=l J 2vr J (/3J) + • • • 



^ E / n [{dA fc }iy [{ A t }]d Wt PK) 



1 

^>0 fc<£ 

J /d0 / A fe (0 / )e^' 7cos( ^^~ a;fe) /-d^" Jd^'A^^Oe^™ 8 ^"-^'"^) 



X £\ / (/3J) •/ 2tt 7 (/3J) 



c 



|{dA}iy[{A}]dcuP(cu) | d0 / A(0 / )e^ Jcos ^-^'-^ + 0(e 2 ) (27) 



27t/ (/3J) 

Thus, with ^(4>) = J{dA}iy[{A}] A(0) we obtain the following (constrained) leading 
order eigenvalue problem, which describes transitions away from the paramagnetic state: 



(0) = ° / ~d<f>' I du P(^)e /3Jcos( ^*'-" ) ^(0 / ) (21 



2vr/ (/?J) 



2tt 







2- 







d0 (0) = (29) 



du; P(w) cos(ku) (30) 



This problem is solved by the Fourier modes ^(4>) = ^e 1 ^ (with integer k ^ 0), each 
of which bifurcates at a temperature which is to be solved from 
1 cI k ((3J) 

WJ) 

For either (3 — > or c — > the right-hand side would reduce to zero, and we would find 
ourselves always in a paramagnetic state. The presently studied transition therefore 
occurs at the average connectivity c for which 

c = min {^7^ F duJ P M cos{ku)\ (31) 



fc>o \I Q ((3J) 

At zero temperature we may use the property linx^oo Ik(z) /Iq(z) = 1 to obtain 
c~ it = max fc>0 fin do; P(c<j) cos(Ax<j) < 1. According to (fT3"|). the present type of 
bifurcation is towards a state where (with denoting the critical Fourier mode): 



1 



i im ITj E W - &]> = IT [! + ecos(A;0 - V) + • • •] (32) 
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Care is to be taken in interpreting the bifurcating state, since in leading order one has 

_e_ / d0 cos(0) cos(A;0 - tp) \ 
2tt ^ jZ n d0 sin(0) cos(A;0 - V) J 

-W28) + - < 33 » 

We conclude that only for k = 1 may we call the bifurcating solution ferromagnetic (F). 
For > 1 we find a bifurcation towards a state with no overall magnetization, but still 
with measurably non-uniform overall single spin statistics. This transition is reminiscent 
of a Kosterlitz-Thouless one. Thus we have the following possible transitions 



J*du P(u) cos(cj) I (34) 



-i 



KT : c = min j ^§ C P(u) cosiku) \ (35) 

It will be shown below, however, that the KT transition is always preceded by a spin- 
glass transition, and hence it is non-physical. 

4-3. Paramagnetic to spin-glass transition 

In those cases where the transition away from the paramagnetic state is towards a new 
state with J{dA}iy [{A}]A(0) = (which would be a spin-glass state), the lowest 
relevant order in our expansions is e 2 , and we find after functional moment expansion: 

/{dA}W[{A}]A(fc)A(fc) = i27r)2 C j2 m f[{dA}W[{A}]du,P(u,)} 

(y d^A^'y^ 08 ^ 1 -*'-^ y d0 , A(0 , )e /3JcosW2 - '-" ) ) + 0(e 3 ) (36) 

Thus, with *(0!,0 2 ) = J{dA}W / [{A}]A(0 1 )A(0 2 ) we now arrive at the following 
(constrained) eigenvalue problem for the P— >SG transition: 

(fc, = cf [2 ^f] )]2 [j^ P( W )^ Jo -^-*i- w H/»Jo-(* a -^- w )j ^ / ) (37) 

J d0x 2 ) = J d0 2 2 ) = (38) 

Again one finds that the Fourier modes are the relevant solutions. Here we have 
*(0i,02) = ^ fcli fc 2 e i(fcl<?il+fc2<fe) (with integer k r ^ 0), each of which bifurcates at a 
temperature (or, equivalently, at a critical connectivity) which is to be solved from 

1 = Clkl ^(pjf J) I dU P{iU) COs[(kl + h)LU] (39) 
As before the right-hand side becomes zero for (3 = or c — > 0, so that the transition 
occurs at 

-.Xl^r/^ PMcosP, + y.]r' (40) 
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Since one always has Ik{/3J) > and I-k(z) = Ik(z), we find that the required extremum 
occurs for k± = — k 2 = k, upon which one subsequently deduces that k — 1, so 

P -> SG : c = ll{(3J)Jll{(3J) (41) 

This equation is obviously independent of the distribution of chiral interactions. At 
zero temperature we find c crit = 1. According to (fTSj) . the present type of bifurcation is 
towards a state where 

1 



Jf E W - hWW ~ &}) = j^y [1+ ecos(0 - 0' + if,) + . . .] (42) 



We note that this particular type of bifurcation obeys liniTv^oo N~ l J2i (fi[4> ~ 4>i\) = 
(2tt)~ 1 , i.e. absence of measurable overall non-uniform spin statistics. Nevertheless, 
limAr^oo N" 1 J2i (^i) 2 > 0, so the bifurcating solution describes a spin-glass state (SG). 

4-4- Summary of transitions away from paramagnetic state and special limits 
We define in accordance with the results (|34l41j) : 

c? 1 = jTdw P(w) cos(c) (43) 

esc = ll(PJ)/Io(PJ) (44) 

The physical transition|| away from the paramagnetic state, as the connectivity is 
increased from c = for fixed /3J, is the one with the largest value of c cr j t . The 
Kosterlitz-Thouless type transition (|3~5j) was given by 



-KT 



max ^ Tin t\ / °^ -^V) cos(A;u;) I < 
fc>i 1 J (/3J) J-* J J (/9J) 



It follows from the properties of the modified Bessel functions |SSj that 4-[l2(z)Io(z) — 
Il{z)\ = \h{z)[h{z) - I (z)} < 0. Hence h(z)/I Q (z) - I 2 (z)/I*(z) < with equality 
only for z — 0. This implies that c<Tq > % T so that the bifurcation (|35)l is indeed 
unphysical. For small and large temperatures one obtains the limiting behaviour 

du P(u) cos(cj) < 1 lim Cp 1 = (45) 

lim Ccp = 1 lim Cop = (46) 

To recover known results in the c — > oo limit one must first re-scale the bond strength 
J. In the case of an overall balance towards (anti-)ferromagnetism, i.e. for chirality 
distributions such that f*^ du P{ui) cos(u;) = O(c ) ^ 0, we have to re-scale according 
to J = J/c. This gives upon taking c — > oo, and using I n (z) = (z/2) n /n\ + 0(z n+1 ): 

T F = ]-J r du P(u) cos(u) (J = J/c, c -> oo) (47) 

Z J — 7T 

Here there is no transition towards an SG state ever. In the absence of such a dominant 
balance, i.e. for distributions such that fl^du P(uS) cos(a;) = A/y/c, the appropriate 

| We will for now leave aside the possibility of first order transitions. 
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c- 1 



SG 



0.0 0.2 0.4 0.6 0.8 1 .0 1 .2 1 .4 1 .6 1.8 2.0 



T/J 

Figure 1. Predicted phase diagram for uniformly distributed chiralities P(uj) = 
(2tt)~ 1 and planar spins (d = 2). It describing a paramagnetic (P) and a spin-glass 
(SG) phase, separated by a continuous phase transition (dotted line). 



re-scaling is J = J '/ ' \fc. Now we find 

T F = l - JA T SG = \j (J= J/V~c, c -> oo) (48) 

Hence as c — > oo and the temperature is reduced from the paramagnetic state, for A > 1 
we will enter a ferromagnetic state, and for A < 1 a spin-glass one. These results are in 
full agreement with those obtained earlier for fully connected systems, see e.g. [30J. 

5. Results for specific chirality distributions at d = 2 

Let us now work out our transition lines (j43|44j) for specific choices for the (symmetric) 
chirality distribution P{uj). It is clear from (|44j) that our choices will only affect the P— >F 
transition. We note that so far we only have expressions for transitions away from the 
paramagnetic state; we are not yet able to determine the F— >SG transition (when both 
F and SG phases exist) analytically, since this would require us to solve our equations 
below the P— >F and/or P— >SG transition temperatures. However we may put forward 
the conjecture (on the basis of our experience with more conventional disordered spin 
models, e.g. that, especially upon taking RSB into account (if needed), there 

will be no change of phase type after the onset of order as the temperature is lowered 
from T = oo to T = 0. This conjecture would predict the elusive F— >SG transition to 
be the horizontal line segment in the (T, c _1 ) phase diagram going from T = to the 
point where the P-^F and P^SG lines meet (the Parisi- Toulouse hypothesis 132J. 

5.1. Predicted phase diagrams 

We first turn to uniformly distributed chiralities: P{oS) = (27r) _1 . Here we find that 
Cp 1 = in (|43|). so for finite c one only ever has the P— >SG transition (|44|). The result 
of numerical evaluation of the latter bifurcation line is shown in figure ^ 
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0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 



T/J T/J T/J 

Figure 2. Continuous phase transitions away from the paramagnetic (P) state for 
planar spins (d = 2) and binary chiralities P(u>) = ^S(ui — ZJ) + \5(uj + uJ), Solid lines: 
P^F bifurcations. Dotted lines: P^SG bifurcations. Note that the location of the 
F^SG transition (dashed) has not been calculated, but follows from the conjecture 
that on lowering temperature the nature of the ordered phase will remain that which 
emerges at the onset. 



A<0 A=\ A=l 




0.0 ' — 1 — 1 — 1 — 1 — ' — 1 — 1 — ' — 1 — ' — 1 — 1 — 1 — 1 — ' — 1 — 1 — ' — 1 — 1 0.0 1 — 1 — 1 — 1 — ' — 1 — ' — 1 — 1 — ' — 1 — 1 — 1 — 1 — ' — 1 — 1 — 1 — 1 — ' — 1 0.0 1 — 1 — 1 — ' — 1 — 1 — 1 — 1 — ' — 1 — ' — 1 — 1 — ' — 1 — 1 — 1 — 1 — ' — 1 — 1 

0.0 0.2 0.4 0.6 0.8 1 .0 1 .2 1 .4 1 .6 1 .8 2.0 0.0 0.2 0.4 0.6 0.8 1.0 1 .2 1 .4 1.6 1.8 2.0 0.0 0.2 0.4 0.6 0.8 1 .0 1 .2 1.4 1 .6 1 .8 2.0 



T/J T/J T/J 

Figure 3. Continuous phase transitions away from the paramagnetic state for planar 
spins (d = 2) and resonant chiralities P{u) = [1 + A cos(£4>)] / 2ir . Solid lines: P— >F 
bifurcations. Dotted lines: P-^SG bifurcations. The location of the F^SG transition 
(dashed) has not been calculated, but follows from the conjecture that on lowering 
temperature the nature of the ordered phase will remain that which emerges at the 
onset. 



Our second choice for the chirality statistics is the binary distribution P(uj) = 
\5{uj — To) + \b~{u) + U), with uJ G [0, 71"]. Here we find equations ()43I44|) reducing to 

_ x _ h(f3J) cos(uJ) _j _ If(/3J) 

Now both types of transition are possible, and it will be clear that this can result in 
richer phase diagrams. The P— >F transitions can only occur if To G [0, 7r/2]; here, as 
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T — > we have < lim^oCF 1 = cos(cu) < lim^oCsG = 1- F° r & = (strictly 
ferromagnetic forces) one will only see a P— >F transition. As we increase ZD away from 
the ferromagnetic value to = we see that at the point where Ii(j3J)/I (pj) = cos(cJ) 
the P— >F transition line crosses the P— >SG one, until at cJ > ^tc the F phase has been 
completely eliminated. The result of numerical evaluation of the bifurcation lines (|4~H|) 
is shown in figure El 

As a third example we inspect resonant chiralities: P{oj) = [l + v4cos(£0)]/27r, with 
A e [-1, 1] and £ G {1, 2, . . .}. Here we find, using Jl n du P(u) cos(cu) = \A5 fA) that 
for A < and also all i > 1 we only have the P^SG transition. If A > and i = 1, on 
the other hand, we may enter either the F or the SG phase: 



2 A WT) Cs " G " WT) m 

with a possible triple point for Ii(/3J)/I (flJ) = \A. The result of numerical evaluation 
of the bifurcation lines is shown in figure El Again we observe the competition between 
ferromagnetic and spin-glass order. 

5.2. Numerical calculation of order parameters via population dynamics 

So far we have only shown results which did not involve solving our order parameter 
equations away from bifurcation points. We now probe our systems further by 
calculating observables (approximately) in the P and SG phases, and by comparing these 
with numerical simulations. Exact execution of this work programme would require us 
to solve the functional W[{P}] from equation (|T6|) . In practice one has to resort to 
explicit parametrizations P[<r|/i] of which the parameters \x are truncated after a finite 
number of components, and solve instead the truncated version of (fTTj) . Here we choose 



m 



exp 


Em>l (^m COS (W>) + ^m sin W)) 




/ d(j)' exp 


E m >i ( A m cos(m0') + A s m sin(m4 


<>))] 



(51) 



with n = (Al, A\, A%, A%, . . .). Using the orthogonality properties of cos(m0) and 
sin(m0) we can extract from ()17j) self-consistent equations for the measure w(fi). For 
the simplest case V(<p) = these equations are found to take the form 

— c I 

w (») = E V" All d» k w(Hk)fakPM] (52) 

£>0 k<£ 

X II d [ A m~2^ l 2n — COs(md>) log e E„>o^, fc ^(^')+^, fc sin(n^)+/3Jcos(0-^. fc ) 
m>l k<l 



Ii^ C m-J2[ — cos(m0) log/ 



x 

m>l k<i 



In order to solve these equations via e.g. the population dynamics scheme [7|, one has 
to truncate the number of coefficients in the parametrization (J51j) . In this paper we 
have limited our analysis to a 2-coefficient truncation, i.e. w(fi) — > w{A\ 1 A\), and to 
populations of size 15. 10 3 in the population dynamics algorithm. Increasing the order of 
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0.5 



q, m 




Figure 4. Population dynamics calculation of observables for uniformly distributed 
chiralities P(u) = (27r) _1 and planar spins (d — 2): q — \{q cc + q ss ) (connected 
open squares) and m = \J m^. + m? s (connected open circles) as functions of c _1 , and 
for T/J = 0.1, 0.3 (upper vs. lower curves). The observables were calculated with 
a population dynamics equation (|52(l and a truncated parametrization (see text for 
details). Connected full squares and circles: corresponding simulation measurements 
of q and m, for systems with N = 10 5 spins. 



the parametrization by small numbers was found to give only a modest improvement in 
accuracy, and it will turn out that already at the present level we find a good agreement 
between theory and (simulation) experiments. 

The numerical simulations, of which data are shown below, were carried out with 
systems of size N = 10 5 and using the Fast Linear Algorithm of [HE]- The measurements 
(values of order parameters) were taken over 10 5 iterations, following an equilibration 
stage of 10 6 iterations. The only exception is FigureEb, obtained using an Euler method 
with elementary time steps of size At = 1/2N and with size N = 1000 (after an 
equilibration stage of 4.10 3 iterations per spin). 

In testing our theory against simulation experiments we focus mainly on the 
following four quantities: 



m. 



N 



N 



N 



Ji^^H ( cos (^)) = / d/x J dcj) P[<p\ii]cos(4>) (53) 

J 1 i ^E( sin (^)) = J^w(^i) |d0 P[0|/x]sin(0) (54) 

i 

Mm ^(cos(0,)) 2 = J dn w(ti) J d0d0' P[0|/x]P[0V] cos(0) cos(0') (55) 

i 

lim = /d/* w(ti) I P[0|/x]P[0V] sm(0) sin(0') (56) 
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Figure 5. Population dynamics calculation of observables for strictly zero chiralities, 
viz. P(ui) = S(ui), and planar spins (d = 2). Left: the scalar observables q = ^(q cc +Qss) 
(connected open squares) and m — \J w? c + ml (connected open circles) as functions 
of c _1 , and for T/J = 0.1,0.3 (upper vs. lower solid curve). Connected full circles 
and squares: simulation data, for N = 10 5 . Right: examples at T = 0.1 of the 
observed distribution P{<p) = N^ 1 ^ ji 5[<p — <j>i\ in simulations (markers), together 
with the corresponding theoretical predictions (solid lines), for c = 10 (left curves) and 
for c = 5 (right curves). All observables were calculated with a population dynamics 
equation 1|52|) and a truncated parametrization (see the main text for details). 



0.8 - 



0.6 - 



q, m 

0.4 



0.2 - 



L 





c- 1 

Figure 6. Population dynamics calculation of observables for binary distributed 
chiralities P(uj) — \5[uj + f ] + \5[uj — j] and planar spins (d = 2): q = \{q cc + q S s) 
(connected open squares) and m = \J m^. + m 2 s (connected open circles) as functions 
of c _1 , and for T/J = 0.2 (upper to lower solid curve). Full markers: simulation data 
(full squares: q; full circles: m), for N = 10 5 . The observables were calculated with 
a population dynamics equation (|52|l and a truncated parametrization (see the main 
text for details). 



These four quantities are then compactified into the following two scalar observables^f: 

<l= -^(Qcc + qss), m = ^m 2 c + m 2 s (57) 

% Many equivalent choices would have been possible. The present definitions have the advantage that 
they will generally give limT^o <l—\ and either limj^o m = (in the SG state) or lim-r^o m = 1 (in 
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These two quantities are sufficient for characterizing any of the three anticipated phases 
{P,F,SG}. In the absence of ferro-magnetism our system is invariant under global 
rotations, which implies that m c = m s = and q cc = q ss = q. Furthermore, in the 
paramagnetic state we have (cos(0j)) = (sin(0j)) = 0, so that we may write 

P : q = m = 
F : q > 0, m > 
SG : q > 0, m = 

The results of our numerical analyses are shown in figures EJ El and |B1 and compared 
with simulation measurements. We find very satisfactory agreement between theory and 
simulation experiments, in spite of the combined limitations imposed by the truncation 
of our parametrization in the population dynamics analysis, the inevitable equilibration 
difficulties of disordered spin systems near their transition points, and the finite system 
sizes in such simulations. 

Figure 13] refers to the chirality distribution P{uj) = (2tt)^ 1 , where we should find 
only the P and SG phases. This is borne out by the data: the order parameter m is 
indeed consistently zero, and q bifurcates to a non-zero value at more or less the predicted 
point. Also the locations of the transitions are as predicted by the corresponding 
phase diagram, i.e Figure In figure El we give data for P(u>) = 5(u). Here all 
interactions are strictly ferromagnetic, leaving only the connectivity disorder, and the 
only possible phases are predicted to be P and F. We do indeed observe the predicted 
non-zero magnetization for small values of c _1 , and again excellent agreement between 
population dynamics and simulations. In this figure we also show the observed and 
predicted shapes" 1 " of the spin angle distribution P((f)) = N^ 1 J2i 8[<f> — (pi] (these are 
predicted by the theory to equal P(4>) = J d/x w(/x)P[0|/x]), for two points in the 
ferromagnetic phase. In the SG phase such measurements tend to be more messy and 
prone to finite size effects, due to the inherent spread of the angles over the interval 
[0, 27t]. As expected, we see that an increase in the connectivity leads to a narrowing of 
the profile of P{4>), reflecting a stronger cooperative ordering of spin orientations. Once 
more the locations of the transitions are in agreement with the phase diagram, i.e. the 
left panel in Fig. 121 Finally, figure El corresponds to the binary chirality distribution 
P{uj) = + f ] + \5[oj — |], and temperature T = 0.2. Here we may test our 

assumption regarding the location of the F— >SG transition line. The prediction of the 
phase diagram (middle panel of Figure EJ) is to have phase F for c -1 < |, phase SG 
for | < c _1 < 0.798133, and phase P for c _1 > 0.798133. The agreement between 
population dynamics and simulations is once more very satisfactory: the magnetization 
m and the spin glass overlap q indeed vanish more or less at the predicted F — > SG and 

the F) state, so that we can always show both quantities together in one plot without loss of clarity. 
+ Since the system is invariant under simultaneous rotations of all spins, the location of the maximum 
of the distribution P(4>) is completely free (only the shape carries information). To enable meaningful 
comparison, one therefore first has to position the theoretical curve such that its maximum coincides 
with that of the simulation data. 
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SG — > P transition points. 

The good agreement between the truncated population dynamics results and our 
numerical simulations underlines the correctness of (i) the RS order parameter equations 
themselves, (ii) the functional moment expansion used to locate the phase transitions 
P^F and P— »-SG, and (hi) the assumed validity in the present models of the Parisi- 
Toulouse hypothesis regarding the location of the F— >SG transition. 



6. General theory for d > 2 

We next try to generalize the theoretical results obtained for d = 2 to d > 2. We will 
define the short-hand \Sd\ — Is^ "- As expected, the paramagnetic solution of equation 
(fTBJl for er G Sd-i and V(tr) = 0, is now seen to be P(<r) = l^-il -1 (as before, it is 
a solution for all T). Our analysis will involve generalizations of the modified Bessel 
functions, such as IqA z ) = (S^-il -1 J s ^dcr e Z(T1 (so that IqA z ) = Iq(z)). 



6.1. Guzai expansion for d > 2 



To find bifurcations away from the paramagnetic state we put P{cr) — > |Sd-i| 1 + A(cr), 
with W[{P}) -> W[{A}] and W[{A}] = as soon as fs^dtr A(<r) ^ so that all 
probability densities are normalized, and we inspect the lowest order functional moments 
J{dA}jy[{A}]A(cr) and J{dA}W^[{A}]A(cr 1 )A(o- 2 ). Close to a continuous transition 
we assume there to be a small parameter e such that / {dAjM/fjA}] A(cr 1 ) . . . A(er r ) = 
0(e r ), as before. We note that 

r /\„/3Jcr.u fe cr'" 



n 

fe=i 



da' 



Sd- 



AJtr') 



\Sd-i\ 



f3.jer-u k er' 



n 



k=l 



1 fs d ^'M<r')e 



J Sd _da'A k (cT , )e> 3J(T - lJ * a " 



1 < / 5d _do-'A fc ( < T')e /3JO '- u '= "/5 d „do- / Afc'(^)e /3J ' T - Ufc ' " 



+ 2^ 

* k^k' 



+ 0(t 



(5f 



Integration of this expression over cr G Sd-i would eliminate the 0(e) term, due to the 
constraint Js^der A(cr) = 0. We may now expand the right-hand side of (|16j) . and find 



£>0 



W[{A}}= Y,-e~ c m^}W[{A k }}dV k P(V k )} J] S 
' ■ k<i cres d _! 



A(o-) + 0{e 3 ) 



i ^ s Sd _d<T'A k (<T>y j <*-w 



\Sd-i 



k<e 



2 \Sd-i\ 1 



ioAPJ) 

f S Sd i d<T'A k {<T'y j(T ^<7' f s ^ i dtr'A k ,(tT')eP Jtr W 



hAPJ) 



hAPJ) 
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\Sd-i\ hAPJ) IoAPJ) 



(59) 



From this follow the relevant functional moment identities. We define *&((t) = 
J{dA}W[{A}}A(cr) and ^(<r\ a 2 ) = J{dA}W[{A}} A^ 1 ) A(<r 2 ). In lowest order e one 
now finds the condition for a continuous P^F transition by solving the (constrained) 
eigenvalue problem 

*(<r) = — <— f ^-M«') I dU P(U)e^™' (60) 

lo,d\pJ ) JS d -i \&d-\\ J 

da *(<t) = (61) 

Sd-i 

If, on the other hand, the first order to bifurcate is e 2 , we find a P^SG transition, 
marked by non-trivial solutions of the (constrained) eigenvalue problem 

^ = W7JT7) I T^*( r ' r2 ) / dU p ^y J((Tl - UTl+<T2 - ur2) (62) 

f da 1 ^(cr 1 , a 2 ) = [ da 2 ^(cr 1 , a 2 ) = (63) 

For d — 2 these expressions ()60I61I62I63|) reduce to those calculated earlier for XY spins, 
as they should. 

6.2. Analysis of the bifurcation conditions 

First we turn to the P^F transition (|60|6ip . It will turn out advantageous to define 
two commuting linear operators K : S^-x — > S^-i and L : S^-i — > S^-i as follows 

(K*)(*) = / -p—eP JtT - T ^{r) (64) 

-'s' d _ 1 |od-i| 

(L*)(tr) = /dU P(U)^(U f cr) (65) 



These definitions allow us to write the eigenvalue problem (|60p61j) as KL^ = 
c ~ l Io,d{fiJ)^ , with constraint f s _ da $?((t) = 0. We can re-write the action of 
L as (LV)((r) = Jda' L(a,a')V(a') with L(a, a') = JdU P(U)5[a' - UV] = 
/ dU P(U)5[a — XJa'}. The kernel L(a, a') then represents the probability that a point 
ct' G Sd-i will be mapped onto a G Sd_i by an orthogonal matrix from the ensemble 
P(U). Both operators K and L are symmetric, as a consequence of -P(U) = P(U'), 
hence we may restrict ourselves to finding simultaneous eigenfunctions of K and L. 

It turns out that a similar strategy can be followed for the P— >SG transition ()62f63|) . 
Here we have to define the commuting linear operators K : Sd ® Sd — > Sd ® Sa and 
L : S d <8> S'd — > Sd <S> Sd as 

(JTO)(<rV 2 ) = / ^^ e ^( - 1 - Tl + CT2 - T2 )*(r 1 ,T 2 ) (66) 
(L\I')(cr 1 , a 2 ) = /dUP(U)*(UV,UV) (67) 
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The eigenvalue problem ()62I63|) can now be written as KL^ = c 1 Jq d (/3J)^, 
with constraints / 5 d<r 1 \l/(er 1 ,er 2 ) = J Sd _ i dcr 2 ^/(cr 1 , cr 2 ) = 0. Here we can re- 
write the action of L as (L^)(cr 1 , cr 2 ) = JdrMr 2 L{o x , cr 2 , r 1 , r 2 )^'(r 1 , t 2 ) with 
Lia 1 , a 2 , t\ t 2 ) = /dU P(U)5[cr 1 -UT 1 ]5[cr 2 -UT 2 ]. The kernel L{cr\ cr 2 , r 1 , r 2 ) now 
represents the probability that a pair of points (t 1 ,! -2 ) G S^-i ® S^-i will be mapped 
onto the pair (cr 1 , cr 2 ) G Sd-i ® by an orthogonal matrix from the ensemble -P(U). 
Both and L are symmetric, as a consequence of -P(U) = P(U^), hence we may once 
more restrict ourselves to finding simultaneous eigenfunctions of K and L individually 

At this stage it would appear appropriate to make a specific choice for the ensemble 
P(U), for which we seek a controlled interpolation between having ferromagnetic and 
completely random chiral interactions. We may define this choice in terms of the above 
probability density L(cr 1 , cr 2 , r 1 , r 2 ) (from which the earlier kernel L(cr, cr') follows by 
integration), for which we now choose the linear combination 

L(cr\ cr 2 , r\ r 2 ) = e^cr 1 - r l ]5[cr 2 - r 2 ] 



+ (1-6) - f J Wn 1 m — i — 91 (6 





-l]5[|cr 2 | 






- i]*[|y| 


- l]<J[a; ■ y - r 1 ■ r 2 ] 



In the non-ferromagnetic part of this measure, i.e. in the term proportional to (1 — e), we 
assign a uniform probability density to all combined image pairs {cr 1 ,^ 2 } of the vectors 
{t 1 , t 2 } which preserve the inner products under the action of the random orthogonal 
matrices U. From definition (|HK|) it then follows automatically upon integration over 
cr 2 that 



L(cr, T ) = eS[cr - r] + (1 - e) - " ' J -. (69) 



5[\cr\ 


-1] 


Is d ^ dx 6 i\ x \ 


-1] 



For e = 1 we return to a strictly ferromagnetic system; for e = we have fully and 
homogeneously distributed random chiral interactions. 

The advantage of our choice (|55j) is that it allows us to diagonalize both kernels 
(166167)) analytically. Working out the eigenvalue problem Js dcr' L(cr, a')^{a') = 
A^(cr) shows that there are only two simple eigenspaces: 

f dcr *(<t) ^ : A = 1 and f dcr ^(cr) = : A = e (70) 

JSd-i J Sd-i 

The first eigenspace is forbidden by constraint (|6ip. so we may simply replace L — > el 
in the eigenvalue problem for the P— >F transition. Working out the P— >SG eigenvalue 
problem J Sd i dr 1 dT 2 L(cr l , cr 2 , r 1 , r 2 )^ (t 1 , t 2 ) = A^(cr 1 ,cr 2 ) leads to the following 

eigenspaces: 

f dcrMcr 2 ^[cr^cr 2 -^^^ 1 ,^ 2 ) = 0, for all w G [—1,1] : A = e (71) 

J Sd-i 

q(cr\cr 2 ) =^(cr 1 - cr 2 ), for all cr 1 , cr 2 G : A = l (72) 

It is clear that the A = 1 eigenspace is perfectly compatible with the constraints (163)) . 
which would e.g. be satisfied by any anti-symmetric function ip{u) in (|72|) . Having 
solved the eigenvalue problem for the operators L for the choice of ensemble (|68)L we 
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may concentrate on the following reduced eigenvalue problems from which to extract 
the phase transitions away from the paramagnetic state: 

P^F: / -^L-e^^ir) = ^M*(«r) (73) 



(74) 



constraint : / dcr ^(cr) = 

J Sd-i 

P SG a : / ^^- e m^ + ^.r^ {T x^ ^ = I lAEl^ l{(T \ ^) 

Js^Sd-xV ec 

f d<r 1 do- 2 «J[<r 1 -<r 2 -M]*(<r 1 ,o- 2 ) = 0, for all u G [—1, 1] 
constraints: f dcr 1 ^(cr 1 , cr 2 ) = / dcr 2 ^(cr 1 , cr 2 ) = 

p _> SGfc . / ^^ e M^.T^^.T 2 )^ (r l. T 2) = 5^)^ ((T 1. ff2) (?5) 

•ASd-i l^d-iT c 
constraint : / dcr ipicr-i) =0 

Although they cannot formally be ruled out, we will henceforth disregard the P^SG a 
transitions, since they are less likely to correspond to the largest eigenvalue in view of 
the extra factor e involved, and because in addition the associated constraints (infinite 
in number) would appear to severely limit the space of allowed functions*. 

6. 3. Explicit results for d = 3 

We finally work out our previous general equations for the value d — 3, where 
we may turn to polar coordinates and write our integration variables as cr = 
(sin(#) cos(0), sin(6 l ) sin(0), cos(6>)) with G [— n, ir] and 9 G [0, 7r]. We note that 
^0,3^) = sinh.(x)/x. Insertion of the polar coordinate representation in our general 
eigenvalue equations shows that for d = 3 the second order P— >F transition is to be 
solved from 

*d6' sin(0') d£ ^(j)^)^)^^)^^^ <p + <p') = sinh ^ J ^ P(g, <p) 



lo 2 J-n 27T \ it -rj ^ J( , 

We observe that the solutions of this equation are of the form \1/(6 I ,0) = p(cos(8))e lk< ^ 
with k G IN (similar to e.g. the spherical harmonics Y^ }jn {6, 0), but with, as we will find 
below, a different dependence on 9), and with the function p to be solved from 

f ^ h{pjjl^yj^)f j *»p{v) = ^T^P(^) (76) 
J-i 2 v epJc 

constraint : either k or J dx p{x) = (77) 

The kernel in (J77|) of which we seek the eigenfunctions is invariant under parity 
transformation, so we may restrict ourselves to eigenfunctions p(x) which are either 

* There are further reasons to reduce the likelihood of the P— >SG a transition being physical. For 
instance, for d — 2 the eigenfunctions can only depend on the inner product between the two vectors 
involved, hence here one simply cannot satisfy the constraints of the P— >SG a bifurcation. 
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symmetric, p+(x), or anti-symmetric, P-{x). Upon implementing these restrictions we 
find 

P+(x) : ['dy I k ((3jVT^JT^) cosh((3 Jxy)p(y) = smh ^ p{x ) ( 78 ) 
Jo v epJc 

p_(x) : fdy I k (/3JVT=tfy/l-y*) smh((3Jxy)p(y) = J) p(x) (79) 

Jo v epJc 

It should be expected that, as for d = 2, the physical (i.e. largest) eigenvalue is the 
one with the lowest allowed value of k, which here is k = 0. Furthermore we note that 
for P-(x) the constraint in ()73|) is automatically satisfied. We have not been able yet 
to solve the above eigenvalue problem analytically, and have instead simply resorted 
to numerical evaluation of the largest eigenvalue in (JBH). This shows that the largest 
eigenvalue is indeed found for k = and P-(x). 

Our analysis of the P— >SG;, transition can be simplified if we use the fact that (|75j) 
is written strictly in terms of inner products of the various unit-length vectors. This 
allows us to choose a convenient basis, e.g. one where cr 2 = (0, 0, 1). Upon again using 
polar coordinates to do the integrations over the sphere S 2 , we find that for d = 3 the 
second order P— >SG transition is to be solved from 

r*d0sin(0) f*d<f> rd6'sm{6') f- d0' eMcosW+xcos(e0+sin(e0cos( ^)^ 



•rr 2,71 Jo 2 27T 

sinh 2 f/? J) 

x V(sin(0) sin(0') cos(0) + cos(0) cos(0 / )) = rn ^t ^) 

c{f3jy 

Via suitable transformations of variables, viz. t = cos(0) and s = cos(0'), and insertion 
of Jdy 5[y — sin(0) sin(0') cos(0) — cos(0) cos(0')] this expression can be simplified to 

1 ^ U U.—W— eM-^ °l(l-s>)(l-t>)-(y-sm 
-i 47r ^f{l-s 2 ){l-t 2 )-(y-st) 2 

sinh 2 (/3J) 

= ^JyT^) (80) 

This latter equation is to be solved subject to the constraint j\dy ip{y) = 0. The 
integration kernel in (J8U|) is again symmetric, yielding as before either strictly symmetric 
or strictly anti-symmetric eigenfunctions. Again the anti-symmetric eigenfunctions offer 
the advantage of automatically satisfying the appropriate constraint. 

We are now in a position to construct phase diagrams for d = 3 and the orthogonal 
random matrix ensemble characterized by 1)68)1 , by solving the remaining two eigenvalue 
problems (J79I80)) with their associated constraints numerically. We have done this for 
three different values of e, leading to the phase diagrams in figure which can be 
compared to the d = 2 results of figure El As might be expected on physical grounds, 
the extra degrees of freedom available to each spin in d = 3 (which increase the potential 
for the system to minimize its free energy entropically, as opposed to energetically) lead 
to a lower transition temperature to an ordered state, be it spin glass or ferromagnetic. 
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T/J T/J T/J 

Figure 7. Phase diagrams in d = 3, for the orthogonal random matrix ensemble 
characterized by equation (|68|l . and for different values of e. e = 1 corresponds to 
purely ferromagnetic interactions, whereas for e = they are fully random. As in the 
case d = 2, except for the location of the triple point, the SG— >F transitions cannot be 
obtained from our present functional moment expansions, but have been inferred from 
the assumption that there is no change of phase character after the onset of order, 
when the temperature is decreased further for fixed connectivity c. 



q, m 




Figure 8. Comparison between theoretical predictions, population dynamics (open 
symbols) and Monte Carlo simulations (solid symbols) for d — 3. Left picture: we show 
the magnetization m (circles) and the spin-glass order parameter q (squares) defined 
in 182(1 . for e = 1 and for connectivity values c = 3,4,5,6 (from left to right). The 
agreement indicates that our truncation in the parametrization of ((81(1 does not have 
a significant impact on the numerical accuracy. Right picture: the order parameters 
m (circles) and q (squares) for e = 1/2, along the line T = 0.25. The magnetization is 
seen to become zero around c _1 = 0.25, thus verifying our assumption for the location 
of the SG-F line. All simulations were done for N = 10 5 spins. 
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6-4- Numerical calculation of order parameters via population dynamics 

Let us now turn to the numerical evaluation of the order parameters in our system. 
In a spirit similar to that of the case d = 2, i.e. section 15. 2| we will extract the 
relevant population dynamics equations from (jl7j) upon making a suitable choice for 
the parametrization of P[<x|/x]. Here we will again use the family of Fourier modes 

P(4> 6\n) = - e XL>i[ A '" cos ( m ^) +A ™ sin ( m ^) + ^ (gl) 

X gSmm'>i[^mm' cos ( m 'l > ) cos { m ' S)+ H^ ml cos (m(f>) sin(m' 9)+ H sin(m</<) cos(m'0)+.ff^ m , sin(m</>) sin(m'S)J 

where D(/j,) is the relevant normalization constant and \i denotes collectively all 
coefficients {A*,B*,H**} with * G {c, s}. From this point, one can proceed 
further by working out expression ()17j) and converting it into one for the coefficients 
{A*, B*, H**}. We have implemented this strategy on the basis of a truncation after 
8 coefficients, i.e. we have taken A* m = B* m = <5 ml and H^ m , = 8 m ^8 m ^\ so that 
H = (A c , A s , B c , B 8 , H cc , H cs , H sc , H ss ) (details of the resulting expressions can be found 
in |Appendix A| ). The self-consistent equation for w(fj,) that follows can then be solved 
using a population dynamics prescription. 

Finally, given the stationary distribution of coefficients u>(/x) in the population 
dynamics method and expression (}8T| one can then evaluate the order parameters of 
the system, e.g. the magnetization and the spin-glass order parameter 

m= ^ m 2 +m 2 + m 2 q= -(q x + q y + q x ) (82) 

where now 

m x = ((cos(0) sin(0))^ |M ) q x = ((cos(0) sin^))^^) (83) 

m y = ((sin(0) siB^))^^) q y = ((sin(0) wn(B))l m ) (84) 

m z = ((cos^))^^} q z = ((cos{6))l m ) (85) 

and with the averages 

((••0W= r# f d0 \sm(6)\ P(J>,9\n) (•••) (86) 

J —TV J 

((.••)> M = /dM^(M) (•••) (87) 

In figure |H1 we show the results of evaluating our observables numerically via the 
above strategy. We compare our bifurcation analysis (leading to predictions for the 
locations of phase transitions) with the outcome of population dynamics analysis and 
with data measured in simulation experiments. In practice, our population dynamics 
shows good convergence already for relatively modest population sizes of 2000 fields. 

The Monte Carlo simulations were done in the present case d = 3 using the 
heat-bath algorithm of (S3 dEj, and with system sizes of 10 5 spins. To generate 
an ensemble of uniformly distributed random orthogonal matrices it is convenient to 
represent the rotation matrices using Euler angles (at, P, 7) (in standard notation), i.e. 
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U = R z (a) Ry(fl) R z {l) f° r an Y U G S0(3). Then, uniform integration over the Lie 
group SO(3) is given by the Haar measure dg#(a,/3, 7) which, in Euler representation, 
takes the form dg#(a, /3, 7) = (87r 2 ) _1 dad/?d7 sin(/3), with a, 7 G [0, 27r] and /3 G [0, 7r]. 
The relevant matrix distribution -P(U) can then be written as 

R z {a) R y {(3) R z {j)) (88) 

(cos 7 sin 7 
— sin 7 cos 7 
1 

In the left picture of figure El we show the magnetization and spin-glass order 
parameters (JH^ for the case where e = 1, i.e. where our rotation matrices reduce to 
the unit matrix U = 1 and the only source of disorder in the system is the random 
nature of the connectivity variables {%}. Here both the population dynamics and the 
simulation experiments reproduce the P-F transition at the predicted connectivity value. 
The excellent agreement obtained indicates that, in retrospect, our truncation of the 
parametrization (}8~Tj) has been made with only a minor cost in accuracy. In the more 
involved scenarios where e < 1, a spin-glass phase with m = and q > will become 
possible. An example is shown in the right picture of figure |Hl where we chose the 
values e = 1/2 and T = 0.25. For these parameter values our theory predicts a P-SG 
transition, but also, based on physical grounds (i.e. absence of re-entrance phenomena) 
we have assumed that the elusive F-SG transition is located at the line segment parallel 
to the T-axis, connecting the triple point where all phases meet to the point T/J=0. In 
figure |H1 we see that also this latter assumption is verified numerically (by the population 
dynamics results), and that the simulation data are, in turn, again in good agreement 
with those of the population dynamics. 

7. Discussion 

In this paper we have applied the equilibrium replica method as developed for finitely 
connected scalar spin systems to models of finitely connected unit-length vectorial spins 
of dimension d, with pair-interactions which are given by random orthogonal d x d 
matrices. Since our spins are continuous and the connectivity c is finite, rather than an 
effective field distribution (as would have been the case for finitely connected discrete 
scalar spins), here the replica-symmetric order parameter is a functional. This generates 
a number of technical complications. Firstly, rather than finding continuous transitions 
away from the paramagnetic state by expansion of the RS order parameter in powers 
of (scalar) moments of a field distribution, here we have to generalize this procedure 
to an expansion in functional moments. Secondly, one should expect serious numerical 
difficulties when attempting to solve numerically the RS order parameter functional 
from the appropriate self-consistent nonlinear population dynamics equation. Here, 
however, these difficulties are kept under control due to the constrained nature of the 



P(U) = e5[U 
with 



1 - e r 2n 



cos a sin a 
R z (a) = —sin a cos a 
1 



/ da dp d 7 sm(/3)S [U - 
Jo Jo Jo 



cos (3 — sin (3 

R y {(3) = \ 1 

sin P cos P 
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microscopic degrees of freedom. Since the present continuous spins live on the sphere 
Sd-x, their value domain is bounded, and it is therefore possible to construct efficient and 
relatively accurate low-dimensional parametrizations of the order parameter functional 
(in contrast to the case of unbounded value domains, as e.g. with ordinary soft spins). 

We have developed our theory initially for arbitrary values of the dimension d of the 
spins, and arbitrary choices of the ensemble of random orthogonal matrices. However, we 
ultimately calculate phase diagrams and the values of moments of the order parameter 
explicitly for d = 2 (where our models reduce to finitely connected XY spins with 
random chiral interactions) and for d = 3 (where they reduce to finitely connected 
classical Heisenberg spins with random chiral interactions). For d = 2, 3 we find three 
types of phases: a paramagnetic (P), a ferromagnetic (F), and a spin-glass phase (SG). 
The calculation of all continuous P— >F and P— >SG transition lines can in all cases be 
reduced to the (numerical) solution of relatively simple functional eigenvalue problems; 
the F— >SG transition is constructed from the location of the triple point, in combination 
with the conjecture (based on previous experience with similar systems, and in line with 
the Parisi- Toulouse hypothesis [32] for the RSB solution) that the phase entered at the 
onset of order (upon leaving the paramagnetic state) will continue to hold upon lowering 
the temperature further for fixed connectivity c. The calculation of observables in the 
F and SG phases was carried out using population dynamics techniques, applied to 
(truncated) parametrizations of the order parameter functional, and the results were 
tested against numerical simulations to reveal excellent agreement. 

We believe the main deliverables of this study to be the successful extension and 
application of finite connectivity replica techniques to more demanding scenarios, where 
the microscopic equilibrated degrees of freedom are neither discrete nor of a scalar 
nature, and where also their interactions are of a mathematically more complicated form 
than just weighted (inner) products. These techniques, which could easily be adapted 
to non-Poissonnian random connectivity graphs, are not only useful tools in the study of 
physical systems, but may also be helpful to analytically determine e.g. the influence of 
topology on global processes in non-physical systems with scale-free connectivity, such 
as the synchronization of randomly and finitely connected planar oscillators |34j . 
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To arrive at a numerically tractable form of the population dynamics relations, one is 
required to truncate ()81|) to a relatively small number of coefficients. Here we took 
1 



exp [A c cos(0) + A s sin(0) + B c cos(#) + B s sin(0)] 



(a.i; 



x exp [H cc cos(0) cos(#) + H cs cos(0) sin(0) + # sc sin(0) cos(#) + sin(0) sin( 

with D(fi) the appropriate normalization constant. The update relations for the above 
coefficients follow from the orthogonality relations of the trigonometric functions: 



ivl 



—c I 

/*) = £ e ~pT / [II d/v(M*)dU fc P(U*)] 



x 5 
x 5 
x 5 
x 5 



irK c (n k ,XJ k ) -4M cs (^ fc ,U A 



7T 7T" 



5 



(A.2) 

nK,^, U fc ) - 4M ss ( ft , U fc N 



k <i <«* - 8) 

7rK j (/x fc ,U fc )-2M(^.,U fc ) 



k<l 



k<l 



7T 7T^ 



F - - £ -M cc (/x fe , U fc )l 4^ CS - E 2 



7rM cs (/x fc ,U fe ) -2# c (/i fe ,U* 



k<t 
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H sc _ -M sc (/z fe , U fc ) 

with the abbreviations 

K c O fe ,U fe ) = 

K s (/x fe ,U fc ) = 

V c (fi k , Ufc) = 

K(Mfc»U fc ) = 

M cc ( Mfe ,Ufc) = 

M cs ( M fc,Ufc) = 

M sc ( Mfc ,Ufc) = 

M ss (/x fe ,Ufc) = 

M(/x fc ,Ufc) = 



k<£ 

E2 

fc<^ 



7r(7T 2 - 8) 

7iM ss (fj, k , Ufc) - 2K s (n k , Ufc) 



7T 



d</> cos(0) / d# 0; Mfc, Ufc) 

7T JO 

"d0 sin(0) fd0 R{(j), 9; fi k , U k ) 

7T JO 

W d<f) [*d6 cos{6) R(<f>, 9; fi k , U k ) 

7T J 

T d0 fd0 sin(0) i?(0,0;/x fc ,Ufc) 

7T JO 



d<j> cos(0) / dfl cos(fl) R(<f>, 9; p k ) 

n JO 

* d<p cosO) 1*69 sin(0) R((f>,6;(i k ,U k ) 

7T JO 

* d0 sin(0) fd0 cos(#) 0; /x fe , U fc ) 

7T JO 

* d0 sin(0) fdfl sin(0) R(<j>, 9; fi k , U k ) 

n Jo 

K d(j) r&e r{4>, e- m„ ^) 

7T JO 



RiM^U) = log(27r)+log 



d9' I sinful e BCco < e '^ BSsia ^ 



(A.3) 



x e 



f3J(U 13 cos(0') sin(e) cos(<f>)+U 2 3 cos(e') sin(0) sin(</>)+£/33 cos(0') cos(0)) 



x J 0, 9', /x, U) + 9, 9', /x, U) 
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Here Iq(x) is the zero-th order modified Bessel function and 

L a (<f>, 6, 6', fi, U) = A c + H cc cos(e') + H cs sm(0') 

+ PJ[U U sin(0) cos(0) sin(0') + U 21 sin(0) sin(0) sin(0') + U 31 cos(6) sin(0') 
L h {<f>, 6, 0', /x, U) = A s + H sc cos(e') + H ss sm(0') 

+ /3j[[/i2 sin(0) cos(0) sin(0') + U 22 sin(0) sin(0) sin(0') + U 32 cos(6) sin(0') 



